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Abstract 


Many  numerical  solvers  aimed  at  the  solution  of  high  speed  chemically  reacting  flows 
employ  methods  that  solve  directly  for  both  the  change  in  thermal  and  mass  diffusivity 
due  to  turbulence  through  the  solution  of  a  turbulent  Prandtl  and  turbulent  Schmidt  number 
directly.  In  many  cases,  the  solution  of  the  turbulent  Schmidt  number  requires  a  method  of 
modelling  the  effects  of  species  massfraction  fluctuations  on  the  source  terms  representing 
the  chemical  reactions  occurring.  One  approach  that  has  seen  reasonable  success  is  the  use 
of  Probability  Density  Functions  (PDFs)  to  evaluate  the  time  averaged  values  of  fluctuating 
massfraction  terms  appearing  in  the  governing  equations.  Both  assumed  forms  for  these 
PDFs  as  well  as  PDFs  where  their  evolution  itself  is  part  of  the  solution  procedure  have 
been  implemented.  The  assumed  forms  have  generally  been  regarded  as  being  the  best 
combination  of  accuracy  and  computational  efficiency.  Indeed,  the  use  of  a  multivariate  /3 
PDF  for  the  species  massfraction  fluctuations  reduces  to  an  algebraic  expression  avoiding 
the  need  to  numerically  integrate  the  PDF  (as  is  required  when  using  a  Gaussian  PDF). 
This  report  details  the  implementation  of  assumed,  jointly  un-correlated  PDFs  for  both 
temperature  and  species  massfraction  fluctuations.  It  shows  that  provided  one  is  already 
using  a  variable  Prandtl  and/or  variable  Schmidt  number  numerical  solver,  no  additional 
conservation  equations  are  required  (the  use  of  a  temperature  PDF  requires  the  variable 
Prandtl  number  solver,  the  use  of  a  composition  PDF  requires  the  variable  Schmidt  number 
solver). 
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Executive  summary 


Report  on  PDF  Models  for  Turbulence  Chemistry 
Interaction 

This  report  describes  the  use  of  Probability  Density  Functions  (PDFs)  in  the  calculation  of 
turbulent,  chemically  reacting  flowfields.  The  requirement  for  the  modelling  of  chemical 
source  terms  is  reviewed  through  the  development  of  the  Favre  averaged  species  conser¬ 
vation  equation.  The  relationship  between  a  time  averaged  and  expected,  or  mean,  value 
(which  can  be  found  using  a  PDF)  is  then  presented.  This  is  followed  by  the  derivation  of 
the  chemical  source  tenn  to  highlight  the  manner  in  which  PDFs  can  be  used  to  solve  for 
this  term  while  considering  fluctuations  in  species  massfractions. 

The  use  of  marginal,  or  jointly  un-correlated,  PDFs  is  presented  using  a  clipped,  Gaussian 
PDF  for  the  fluctuations  in  temperature  and  a  multivariate  /3  PDF  for  species  massfraction 
fluctuations.  It  is  shown  how  the  use  of  the  multivariate  /3  PDF  reduces  to  the  solution  of 
a  set  of  algebraic  equations  eliminating  the  need  to  numerically  integrate  the  PDF.  This 
algebraic  expression  requires  the  solution  of  the  time  averaged  species  massfraction  and 
the  sum  of  the  species  massfraction  variance,  where  the  variance  is  assumed  known  (this  is 
often  solved  for  as  part  of  a  variable  Schmidt  number  numerical  solver). 

The  report  also  contains  a  brief  review  of  relevant  papers  utilizing  the  methods  described 
in  this  work  and  concludes  with  some  recommendations  for  future  work. 
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1  Introduction 


The  determination  of  fluctuating  components  within  a  turbulent  flowfield  has  long  been 
an  area  of  significant  computational  research  given  the  ease  with  which  laminar  values,  or 
time  averaged  values,  can  be  calculated  using  reasonably  efficient  Navier-Stokes  solvers. 
This  has  led  to  a  wide  variety  of  turbulence  models  whose  aim  is  to  accurately  model  the 
effects  of  fluctuations  in  the  velocity  field  on  the  apparent  viscosity  within  the  flowfield. 
However,  velocity  is  not  the  only  variable  to  fluctuate  due  to  turbulence  and  thus  numerous 
researchers  have  led  efforts  into  the  direct  modelling  of  other  fluctuating  flowfield  variables. 

Work  done  by  Nagano  and  Kim  [  1  ]  introduces  a  two  equation  turbulent  thermal  field  model. 
Sommer,  So,  and  Lai  [2],  So  and  Sommer  [3],  and  Sommer,  So,  and  Zhang  [4], [5]  build 
on  this  work  by  modifying  the  turbulent  thermal  field  equations  to  include  the  effects  of 
viscosity  and  compressibility  while  modifying  the  near  wall  treatment  to  obtain  behaviour 
consistent  with  the  physics  of  the  flow.  This  results  in  a  numerical  method  that  is  able  to 
calculate  the  turbulent  Prandtl  number  directly  without  a  need  to  set  it  a  priori,  allowing  a 
more  accurate  simulation  of  flows  in  which  this  number  can  vary  significantly  within  the 
flowfield  (like  rocket  plumes  or  scramjet  combustors). 

For  multi-species  flows  turbulence  can  increase  the  apparent  mass  diffusivity  (just  as  it  does 
for  both  viscosity  and  thennal  diffusivity)  which  has  led  to  the  development  of  variable 
Schmidt  number  numerical  solvers.  Brinkman  et  al.  [6]  consider  non-reacting  flows  and 
propose  that  the  turbulent  Schmidt  number  can  be  calculated  directly  from  the  solution 
of  an  equation  governing  the  sum  of  the  square  of  the  fluctuations  in  the  species  mass 
fractions  and  another  governing  the  rate  of  this  quantities’  dissipation.  This  approach  is 
further  expanded  by  Brinkman  et  al.  [7]  [8]  and  Calhoon  et  al.  [9]  where  the  consideration 
of  reacting  flows  is  included  (while  also  replacing  the  thermal  variance  equation  with  one 
for  the  internal  energy  variance).  In  each  of  these  cases  these  variable  Schmidt  number 
models  are  coupled  to  a  two  equation  k  —  (0  turbulence  model.  It  is  interesting  to  note  that 
in  these  works  even  though  chemical  reactions  are  considered,  the  need  to  model  chemical 
reaction  source  terms  is  avoided  through  the  use  of  two  conservation  equations,  one  for  the 
mixture  fraction  variance  and  another  for  the  rate  of  its  dissipation. 

Although  not  as  common  as  the  two  equation  k  —  (0  turbulence  model,  there  has  also  been 
development  of  a  variable  Schmidt  number  model  with  the  k—  £  turbulence  model  (Robin¬ 
son  and  Hassan  [10]).  This  model  has  been  found  to  work  well  when  using  both  a  variable 
turbulent  Prandtl  and  Schmidt  number  approach  (Xiao  et  al.  [1 1][12][13]). 

In  general,  when  solving  for  chemically  reacting  flows  the  use  of  a  species  massfraction 
conservation  equation  leads  to  the  need  to  model  the  average  of  a  chemical  source  term. 
Terms  involving  the  products  of  species  massfraction  fluctuations  and  the  chemical  source 
term  appear  in  other  equations  as  well  and  thus  require  additional  modelling.  A  popular 
approach  in  handling  these  tenns  has  been  through  the  use  of  probability  density  functions, 
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or  PDFs.  Two  popular  approaches  for  this  method  have  been  to  use  either  an  assumed 
form  for  the  PDF,  or  to  calculate  the  evolution  of  the  PDF  as  part  of  the  solution.  The 
evolution  approach  is  based  of  the  work  of  Pope  [14],  [15]  where  the  evolution  of  a  joint 
velocity/composition  PDF  is  solved  using  a  Monte  Carlo  approach.  The  advantage  of  this 
method  is  that  the  effects  of  convection,  reaction,  body  forces,  and  the  mean  pressure  gra¬ 
dient  appear  directly  in  the  equations  and  thus  do  not  require  modelling.  However,  due  to 
the  large  dimensionality  of  the  joint  PDF  this  method  is  computationally  expensive  and 
thus  many  researchers  have  focused  on  using  an  assumed  form  for  the  PDF  based  on  the 
work  of  Girimaji  [16],  [17]  called  the  multivariate  /3  PDF. 

When  comparing  the  two  methods  both  have  been  found  to  produce  similar  results  for  the 
mean  flow  variables,  however,  the  evolution  PDF  methods  require  a  significant  increase  in 
computational  requirements  (Calhoon  and  Kenzakowski  [18],  Baurle  et  al.  [19]  [20]  [21]). 
In  all  these  cases  the  temperature  fluctuations  are  modelled  using  an  assumed  Gaussian 
distribution  (or  assumed  /3  function  for  Calhoon  and  Kenzakowski)  while  the  massfraction 
fluctuations  are  modelled  using  the  assumed  multivariate  /3  PDF  of  Girimaji.  The  assumed 
multivariate  /3  PDF  approach  has  difficulty  calculating  higher  order  terms  due  to  its  as¬ 
sumption  of  statistical  independence  between  the  temperature  and  the  composition.  It  is 
also  fairly  dissipative  in  the  presence  of  combustion  (Keistler  et  al.  [22])  making  it  less 
accurate  downstream  of  ignition,  but  reasonably  accurate  preceding  ignition.  It  has  been 
shown  to  be  fairly  insensitive  to  the  choice  of  turbulence  model  in  terms  of  its  effect  on 
the  calculation  of  the  turbulent  mass  diffusivity  and  requires  no  additional  variables  to  be 
solved.  The  requirement  of  no  additional  unknown  terms  other  than  those  already  present  in 
most  variable  Schmidt  number  solvers  makes  the  assumed  multivariate  /3  PDF  an  attractive 
choice  for  many  computational  codes.  Due  to  the  temperature/composition  de-coupling  it 
is  also  possible  to  treat  the  fluctuations  in  temperature  and  composition  differently. 

Narayan  and  Girimaji  [23]  use  an  assumed  form  multivariate  /3  PDF  to  model  massfraction 
fluctuations  while  avoiding  a  PDF  for  temperature  through  the  use  of  a  series  expansion  for 
the  reaction  rates  based  on  the  mean  temperature.  Alternatively,  Gaffney  et  al.  [24]  ignore 
the  effect  of  massfraction  fluctuations  and  model  temperature  fluctuations  as  following  ei¬ 
ther  an  assumed  Gaussian  or  assumed  j3  PDF.  They  show  that  temperature  fluctuations  can 
have  a  significant  impact  on  the  calculation  of  reaction  rates,  while  the  use  of  a  Gaussian 
PDF  can  be  sensitive  to  the  clipping  limits  used  (since  the  exact  Gaussian  PDF  has  limits 
of  ±°°).  In  addition,  Gaffney  et  al.  [25]  have  studied  series  expansion  for  the  fluctuations 
in  composition  with  assumed  PDFs  for  temperature  resulting  in  a  method  that  does  not 
require  any  numerical  integration.  Tested  on  a  high  speed  turbulent  reacting  hydrogen/air 
mixing  layer  (two  dimensional)  they  show  that  species  fluctuations  can  also  have  a  notice¬ 
able  effect  on  the  chemical  reaction  rates. 

The  multivariate  /3  distribution  of  Girimaji  for  the  species  massfraction  PDF  has  also  been 
tested  with  various  assumed  temperature  PDFs  by  Gerlinger  [26].  The  shape  of  the  assumed 
temperature  PDF  (Gaussian,  triangular,  rectangular)  is  found  to  have  little  effect  on  the  re- 
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suits  of  a  supersonic  diffusion  flame.  However,  the  choice  of  variance  equation  from  which 
the  temperature  fluctuations  are  extracted  shows  that  a  semi-empirical  equation  achieves 
the  best  results  (where  transport  equations  for  the  variance  of  sensible  energy,  energy,  and 
temperature  are  compared). 
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2  Turbulent,  Reacting  Flow 


For  a  system  in  which  the  fluid  composition  can  change  due  to  chemical  reactions,  it  is 
often  convenient  to  model  a  conservation  equation  for  each  of  the  possible  species, 


dt 


(pYm)  "F  (pw/fm) 


dxj 


pD„ 


djrn 

'  dxj 


dt;i;  —  0 


where  Ym  is  the  massfraction  of  species  m, 


(1) 


Ym  —  Pm/ P  (2) 

Although  accurate  for  laminar  flows,  Eq.  1  does  not  account  for  the  fact  that  under  turbu¬ 
lent  conditions  each  of  the  variables  may  fluctuate  at  a  given  point  and  thus  under  these 
circumstances  an  average  over  time  can  be  taken, 


dt 


(pYm)  4-  ( pUjYn/j 


dXj 


pDn 


dYm\  _- 


dXj 


(dm  —  0 


(3) 


where  the  averaging  process  indicated  by  the  overline  is  often  referred  to  as  the  Reynolds 
average.  For  a  given  variable  (f), 


—  1  ft=T 

<P=r  Qdt  (4) 

1  Jt= o 

It  is  also  possible  to  define  a  density  weighted  average  (often  referred  to  as  a  Favre  average) 
as, 


$ 


1 

w 


p± 

p 


(5) 


where  it  is  possible  to  relate  several  Reynolds  and  Favre  averaged  quantities  through, 


P<P  =  P<t> 


(6) 


and 


P0102  =p0l02  +  P<K>2 


(7) 
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Therefore,  if  one  assumes  a  variable  can  be  represented  by  its  Favre  averaged  value  and  a 
fluctuation  about  this  value,  then  using  species  massfraction  as  an  example, 


and  Eq.  3  can  be  expressed  as, 


(8) 


In  obtaining  this  expression  it  has  been  assumed  that  the  dominant  effect  of  turbulence 
on  the  mass  transfer  process  is  from  the  large  scale  movement  of  turbulent  eddies  (thus 
allowing  one  to  assume  the  effects  of  u"  to  be  much  greater  than  those  due  to  D" ).  In 
this  form  the  solution  becomes  one  of  finding  the  time  averaged  values  of  quantities  such 
as  velocity,  species  massfraction,  and  density.  However,  this  inclusion  of  turbulence  now 
involves  the  solution  of  an  additional  term  involving  the  fluctuating  components  of  both 
velocity  and  massfraction,  as  well  as  a  time  averaged  value  of  the  species  production  6)m. 

Most  numerical  models  employ  the  assumption  of  a  turbulent  diffusivity,  ( vm )  r  =  ~p(Dm)T, 
which  is  used  to  approximate  the  fluctuating  terms  through, 


(9) 


allowing  the  species  conservation  equation  to  take  the  form, 


The  calculation  of  the  turbulent  diffusivity  is  often  based  off  of  a  knowledge  of  the  turbulent 
co-efficient  of  viscosity  jj.r-  This  term  accounts  for  the  apparent  increase  in  viscous  effects 
due  to  the  fluctuations  in  velocity  and  is  used  to  model  the  Reynolds  stresses  which  appear 
in  the  momentum  equations  after  performing  a  similar  averaging  process,  i.e., 


(11) 


There  are  numerous  approaches  to  calculating  a  value  for  fir-  many  of  which  involve  solv¬ 
ing  an  additional  one  or  two  conservation  equations  for  various  turbulent  quantities.  Once 
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found,  if  one  sets  a  turbulent  Schmidt  number  a  priori  then  this  leads  directly  to  a  value  for 

(Dm)T  as, 


(Dm)T  =  ^~  (12) 

pScT 

Although  this  approach  can  produce  satisfactory  results,  it  is  dependent  on  a  knowledge  of 
the  correct  turbulent  Schmidt  number  to  set  for  the  particular  flowfield  under  consideration 
(and  thus  must  often  be  tailored  for  a  given  simulation).  However,  it  is  also  possible  to 
directly  calculate  the  turbulent  diffusivity  in  a  manner  similar  to  that  done  for  the  turbu¬ 
lent  co-efficient  of  viscosity  through  the  solution  of  additional  conservation  equations.  In 
either  case,  the  sole  remaining  task  for  the  solution  of  Eq.  10  is  some  means  of  finding  the 
averaged  value  of  the  species  production  tenns. 

2.1  Expected  Value 

In  the  previous  section  the  assumption  of  a  time  average  is  used  to  postulate  the  existence 
of  a  value  that,  given  a  large  enough  time,  will  be  consistently  measured  for  a  given  set 
of  circumstances  (i.e.,  initial  conditions).  This  implies  that  despite  the  random  nature  of 
turbulence,  the  averaged  values  solved  for  would  not  vary  significantly  between  successive 
repetitions  of  the  same  experiment.  The  variability  across  successive  experiments  (each 
conducted  over  a  given  time  period)  is  evaluated  by  the  ensemble  average, 

(13) 

iV  N 

where  the  number  of  successive  repetitions  of  the  experiment  is  N.  As  the  number  of  rep¬ 
etitions  increases,  the  value  of  the  ensemble  average  tends  towards  the  mean,  (</>),  or  the 
expected  value,  which  is  formally  defined  as, 


/oo 

y/P^(y/)dy/  (14) 

-oo 

In  Eq.  14  PA y/)  is  called  the  probability  density  function  and  it  represents  the  probability 
that  the  variable  (f>  will  have  a  value  equal  to  y/,  where  y/  represents  all  the  values  possible 
for  the  variable  (p  (i.e.,  all  the  realizations  of  (p).  In  the  case  of  a  Gaussian  random  variable 
(also  known  as  a  normally  distributed  random  variable),  the  probability  density  function 
(hereafter  referred  to  as  the  PDF)  has  a  bell  shape  defined  by, 


p*(v) 


y/2no 


(15) 
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where  o2  is  the  standard  deviation.  Depending  of  the  values  of  the  mean  and  the  standard 
deviation  the  resulting  curve  can  be  shifted,  sharpened/flattened,  or  both  (see  Fig.  1) 


Figure  1:  Gaussian  probability  distribution  function 


If  the  exact  same  PDF  applies  across  multiple  experiments  (i.e.,  not  simply  the  shape  of 
the  PDF  but  also  the  mean  and  standard  deviation)  then  the  experiments  are  said  to  be 
statistically  stationary  (they  have  the  same  probability  law)  and  thus  from  the  definitions  in 
Eq.  13  and  14, 


(0)tf  —  jy  K0)i  +  (0) 2  H - F  (0)w] 


^n=N 


yP^(y)d  y  +  /  yP<j>2(y)dy+---  +  /  yP(j)N(y)dy 


but  since  p*(v)=p*  (y)  =  P$N(y)  this  becomes, 


(<I>)n 


1 

N 


yP<i)(y)dy 


=  (0) 


(16) 


for  all  values  of  N,  including  N  —  \  .  Thus  for  statistically  stationary  flows  the  ensemble  av¬ 
erage  is  exactly  equal  to  the  mathematical  expectation  even  for  a  single  set  of  experimental 
results. 


If  one  postulates  that  as  the  time  interval  over  which  an  average  is  taken  increases  (larger 
T),  the  difference  between  values  obtained  (as  defined  by  Eqs.  4  or  5)  will  become  smaller 
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between  successive  repetitions  of  a  given  experiment,  then  the  time  averaged  value  (</)) 
approaches  the  ensemble  average  ((0)#).  If  the  flow  is  further  assumed  statistically  sta¬ 
tionary  then  by  virtue  of  Eq.  16  this  means  that  the  time  averaged  value  of  a  given  variable 
approaches  the  mathematical  expectation,  or  mean,  of  that  variable, 


0  =  {«  (17) 

This  result  provides  an  interesting  approach  to  finding  the  time  averaged  value  of  any  vari¬ 
able  (including  the  time  average  of  a  fluctuating  variable)  through  the  use  of  PDFs.  Before 
examining  this  approach  in  detail,  it  is  convenient  to  define  the  nth  central  moment  (or  the 
nth  moment  about  the  mean)  as  the  expectation, 


/oo 

ty -{$)]"  P${v)dv  (18) 

-oo 

where  since  the  mean  lies  at  the  centre  of  the  PDF  the  first  central  moment  is  equal  to  zero. 
However,  the  2nd  central  moment  (also  referred  to  as  the  variance)  is  directly  related  to  the 
standard  deviation,  G^,  as  the  expectation  of  the  square  of  obtaining  a  value  away  from  the 
mean, 


var{(j))  =  oj=E  (ly/-^)]2^  =Jj¥-m2P$Md¥  (19) 

This  provides  a  link  between  fluctuations  about  a  mean  value  and  the  PDF.  For  any  variable 
composed  of  a  time  averaged  and  fluctuating  component,  when  statistical  stationary  flow 
is  assumed  one  can  write, 


(j)  =  $  +  f  =  ((j))  +  f  =>$"  =  $-  (<p) 

The  expected  value  of  this  fluctuating  quantity  squared  can  thus  be  expressed  as, 

£(O")=£([0-<*>]2)  (20) 

where  comparing  the  righthand  side  of  Eq.  20  to  Eq.  19  shows  that  the  standard  deviation 
is  directly  related  to  fluctuating  component  of  the  variable  under  consideration, 

E  (O")  =  <0">  =  a}  (2D 
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2.2  Chemical  Source  Term 


The  solution  to  Eq.  10  requires  a  value  for  6)m,  the  time  averaged  rate  of  production  or 
consumption  of  species  m  per  unit  time  ([(kg/(m3/s)]).  For  laminar  flows  this  rate  can  be 
calculated  from  a  knowledge  of  the  forward  and  backwards  reaction  rate  of  a  given  chem¬ 
ical  reaction  along  with  the  stoichiometric  mol  numbers  of  a  given  species  within  that 
reaction  as  follows.  Using  the  dissociation  of  hydrogen  through  a  collision  with  another 
molecule  ( B )  as  an  example, 


one  can  define 


2  H  +  B^H2+B 


(22) 


and 


vr  ■ 

..reactant  0 

vH  — z 

..reactant  A 
vHl  -0 

..reactant 

VB 

yP  ■ 
ym  * 

^product  =  Q 

product  9 
vh2  —  1 

^product 

(23) 


For  a  molecule  that  does  not  undergo  a  change  in  quantity  through  the  reaction  the  differ¬ 
ence  between  —  v/„  will  always  be  zero  while  for  one  which  increases  in  the  forwards 
direction  this  difference  will  yield  a  positive  value  (and  thus  a  negative  value  is  obtained  for 
a  molecule  that  decreases  in  the  forwards  direction).  Therefore,  if  one  assumes  an  overall 
reaction  rate  K  such  that  a  positive  value  indicates  the  forward  direction  (from  left  to  right 
in  Eq.  22)  then  the  chemical  source  term  can  be  expressed  as, 


d>m  =  ( <  -  v’m  )K  (24) 

provided  the  overall  reaction  rate  K  has  units  of  (m3s)_1 .  This  rate  is  found  as  the  difference 
between  the  rate  at  which  the  given  reaction  proceeds  in  the  forwards  and  backwards.  The 
forward  reaction  rate  can  be  calculated  using  a  modified  Arrhenius  equation, 

kf=ATne~Ea/^T  (25) 

and  the  backwards  reaction  rate  is  related  to  kf  through  the  equilibrium  constant  based  on 
concentrations, 


(26) 


It  should  be  noted  here  that  by  definition  the  equilibrium  constant  based  on  concentrations 
is  defined  as, 
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KC(T) 


)v'" 


n  m 
n  m 


(pYm\ 

\  -  YYm  J 


( pYm\ 

\  -  Yf  >n  J 


Kc(T)  =  n,„ 


(27) 


where  from  Eq.  27  it  can  be  seen  that  the  equilibrium  constant  will  have  units  of  concen¬ 
tration  (kmols/m3)  raised  to  the  power  of  X  Y  —  X  v/„ .  This  requires  that  the  various  terms 
appearing  in  Eq.  25  (activation  energy  Ea  (per  kmol),  the  exponent  n,  and  the  constant  A) 
which  vary  depending  on  the  particular  molecular  reaction  being  considered  (91  is  simply 
the  universal  gas  constant)  be  such  that  they  yield  a  rate  (i.e.  change  per  unit  time)  divided 
by  concentration  raised  to  the  power  of  the  sum  of  the  reactant  mol  numbers,  Xm  v'r  From 
Eq.  26  this  will  yield  a  backwards  reaction  rate  with  units  of  per  unit  time  and  per  unit 
concentration  raised  to  the  power  of  the  sum  of  the  product  mol  numbers  (Xvm)-  Under 
these  circumstances,  the  overall  reaction  rate  can  then  be  calculated  on  a  simple  per  unit 
time  basis  by  a  multiplication  of  the  forward  rate  by  the  sum  of  the  reactant  concentrations 
(each  raised  to  their  stoichiometric  mol  numbers,  see  Eq.  23)  minus  the  backwards  reaction 
rate  multiplied  by  the  product  concentrations  (raised  to  their  stoichiometric  mol  numbers). 


K  —  k 


f 


pYHy»  / pY^Y*  (pYB 

.-Mj-j  J  \  'h2  J  V 


-kb 


pyh\v»  / pYH2y»2  (pYb 

■  J  \^H2) 


=  i-,nfr&+X&,+XS) 


K  =  kfp 
~hp 


Yh_YH  ( I?h_YH2  flB 

#H2 


<+ypH0+<)  (  Ih_YPh  ( Ith_Y»2  (jn\v^ 

%h2 


Defining  the  variables 


C  =  Iv;  and  {'’  =  £' 


(28) 


allows  one  to  re-express  the  chemical  source  tenn  in  Eq.  24  as, 


(Dm  =  -JZmiVm  -  v„) 


kfp?  nr 


v 

1  m 


-kbp^’nn 


Y 

1m 
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In  the  case  where  there  is  more  than  one  reaction  that  involves  the  species  m,  the  above 
expression  is  modified  slightly  (where  the  number  of  possible  reactions  ranges  from  1  to 

J)> 


CO, 


bm  -y$rn  ' 


'vp  --V  ■] 

<  mj  vm,j J 


k/.pV  n 


n, 


y\  ’W2,/ 
m  x 


(29) 


as  now  the  stoichiometric  mol  numbers  of  species  as  products  and  reactants  depends  on  the 
particular  reaction  j  and  thus, 


(j=Kj  and  ff  =  I  <j  (30) 

m  m 

while  the  definitions  of  kf  and  kb  are  the  same  as  those  in  Eqs.  25  and  26  (only  now  specific 
to  reaction  j). 

With  Eq.  29  one  is  now  in  a  position  to  take  the  time  average  of  the  chemical  source  tenn, 
where  since  the  molecular  weights  and  the  stoichiometric  mol  numbers  are  constant  results 
in  the  following, 


COiij  — 


m,j 


~  V, 


mjJ 


kfjp^'  n m  ( *^7 


- V 

nJY  m’J 

J  1  m 


hjp  n  m 


— v 


I:j  Y  rnj 
J  J-  m 


In  this  form  the  challenge  becomes  how  to  properly  evaluate  the  average  reaction  rates 
(both  forward  and  backward)  along  with  the  fluctuations  in  species  massfraction.  For  ex¬ 
ample,  including  turbulent  fluctuations  in  temperature  the  forward  reaction  rate  (Eq.  25) 
yields, 


kfj=A(T+T")ne^T+T")  (31) 

where  even  if  T" ft  is  assumed  <  1  thus  allowing  (T  +  T")  to  be  expanded  in  a  series  still 
requires  numerous  terms  to  be  carried  through  the  calculations  to  be  accurate.  Similarly, 
when  including  turbulent  fluctuations  in  species  massfraction  one  obtains  the  terms, 


n m{Ym  +  Y£)  m'J  =  (Y1+Y;/)^J(Y2  +  Y^J---(Ym  +  Y^j 


which  results  in  a  large  number  of  fluctuating  massfraction  products  that  would  require 
modelling. 


11 


Therefore,  in  order  to  avoid  these  difficulties  one  can  replace  the  average  value  for  kj.  kb, 
and  Ym  directly  with  their  mean  values  if  one  assumes  statistically  stationary  flow  (Eq.  17), 


6)m  — 

(32) 

It  should  be  noted  that  in  obtaining  Eq.  32  it  has  been  assumed  that  the  density  and  the 
massfractions  are  jointly  independent  (see  Gaffney  et  al.  [27]). 

2.3  Joint  PDFs 

When  using  PDFs  to  find  the  mean  value  of  a  variable,  if  the  variable  is  found  to  depend 
on  more  than  a  single  independent  quantity  then  a  joint  PDF  is  required.  However,  this 
then  requires  that  to  evaluate  the  joint  PDF  one  must  be  capable  of  integrating  over  all  the 
independent  quantities.  For  the  case  of  the  chemical  source  tenn,  the  rate  at  which  a  given 
species  is  produced  or  consumed  depends  on  both  the  temperature  of  the  mixture  as  well 
as  the  massfractions  of  the  various  species  that  exist  at  a  given  time.  This  means  that  a  joint 
PDF  of  both  temperature  and  composition  is  required, 


(33) 

As  with  a  single  variable  PDF,  when  the  integration  of  the  joint  PDF  is  taken  over  the 
limits  of  negative  and  positive  infinity  the  results  is  still  unity  (i.e.,  considering  the  entire 
spectrum  of  possible  values  for  a  given  variable,  the  probability  is  certain), 


/  ■■■/  .,&m)(T,Y\, ...  ,Ym)dTdYi  ■■■dYm  —  1  (34) 

However,  if  one  assumes  that  the  species  massfractions  are  jointly  independent  of  the  tem¬ 
perature  then  the  joint  PDF  can  be  expressed  as  a  multiplication  of  two  separate  PDFs, 


u...,&m){T iY\  P^(T)P^U^(YU. . . ,  Ym) 


each  of  which  can  be  evaluated  independently. 


(35) 
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2.4  Assumed  Temperature  PDF 

Although  the  assumption  that  the  temperature  and  species  massfractions  are  uncorrelated 
(jointly  independent)  does  not  reflect  the  physics  of  turbulent  combustion,  it  has  been  ap¬ 
plied  successfully  to  mixing  flows  and  thus  extended  to  flows  involving  combustion  as  well. 
When  this  is  done,  the  temperature  can  be  assumed  to  have  a  randomly  distributed  profile 
and  thus  follow  a  Gaussian  PDF, 


PAT) 


y/2nor 


(36) 


For  statistically  stationary  flows  one  can  replace  (T)  with  T  while  from  Eq.  21  the  standard 
deviation  of  the  PDF  can  be  replaced  with  the  average  of  the  square  of  the  fluctuations  in 
temperature, 


]  ( T-T )2 

Pg-{T)  =  -  =e  2 t"t" 

V2kT"T" 


(37) 


Although  the  Favre  averaged  temperature  is  likely  a  known  variable  through  the  standard 
solution  of  the  Favre  averaged  energy  equation,  the  fluctuating  temperature  is  generally 
not  a  variable  in  the  solution  vector.  However,  several  authors  have  developed  methods  to 
avoid  the  reliance  on  a  pre-determined  value  for  the  turbulent  Prandtl  number  allowing  the 
direct  calculation  of  the  turbulent  eddy  thermal  conductivity  ( Kj )  from  a  set  of  equations 
similar  to  those  used  for  solving  for  the  turbulent  eddy  viscosity  (/J-t).  These  methods,  often 
referred  to  as  variable  Prandtl  number  models,  generally  solve  for  a  variance  equation  (i.e., 
the  expectation  of  the  square  of  a  fluctuating  value)  and  its  dissipation.  For  use  in  evaluating 
the  turbulent  eddy  thermal  conductivity  the  variance  of  either  temperature  (Nagano  and 
Kim  [1]),  specific  enthalpy  (Xiao  et  al.  [28], [13]),  or  specific  energy  (Brinckman  et  al.  [7]) 
have  been  modelled.  In  cases  where  the  temperature  variance  is  not  modelled  directly  it  can 
be  recovered  by  neglecting  the  effect  of  turbulent  fluctuations  on  the  specific  heat  allowing 
one  to  write, 


T"T"  =  (^CPmYm)  h"h"  or  T" T"  =  e"e"  (38) 

\  m  J  \  m  J 

Therefore,  most  variable  Prandtl  number  approaches  will  already  include  sufficient  infor¬ 
mation  to  completely  define  the  temperature  PDF  in  Eq.  37  and  thus  the  mean  value  of  a 
given  reaction  rate  (either  forward  or  backward)  can  be  found  through, 


poo 

(kfj)  JQ  kfjP^(T)dT 


(39) 
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Several  items  should  be  noted  when  using  Eq.  39.  The  first  is  that  the  limits  on  the  integra¬ 
tion  no  longer  extend  below  zero.  Physically  this  is  justified  in  that  the  temperature  cannot 
fall  below  this  value,  however,  a  truly  Gaussian  PDF  requires  that  whatever  value  is  being 
considered  has  no  limits  so  that  the  requirement  in  Eq.  34  is  respected, 


Psr{T)dT=  1 


(40) 


Furthermore,  when  numerically  integrating  Eq.  37  finite  limits  on  both  ends  must  be  ap¬ 
plied  leading  to  a  further  clipping  of  the  Gaussian  PDF.  These  limits  should  be  set  so  as  to 
obtain,  as  close  as  possible,  the  result  in  Eq.  40  over  the  range, 


P^(T)dT^  1 


(41) 


while  also  respecting  the  fact  that  over  this  range  the  forward  reaction  rate  must  be  cal¬ 
culated  according  to  the  model  constants  in  Eq.  25  which  are  valid  within  a  finite  range 
(generally  between  ~  300  and  ~  3000  K  depending  on  the  model  chosen).  Thus  the  final 
result  for  the  mean  reaction  rate  can  be  expressed  as, 


kf,(T)  - 

-  - P 

V2kT77T77 


( T-f )2 

2  T"T"  dT 


(42) 


Of  course,  if  the  variance  is  zero  then  the  Gaussian  PDF  approaches  the  behaviour  of  a 
delta  function  centred  at  the  mean  which  can  cause  difficulties  for  numerical  integration. 
However,  as  the  variance  approaches  zero  this  means  that  the  turbulent  fluctuations  are 
approaching  zero  and  thus  the  flow  can  be  treated  as  laminar,  in  which  case  the  reaction 
rate  calculated  directly  from  Eq.  25  can  be  used  with  T  —  T . 

2.5  Assumed  Massfraction  PDF 

To  determine  the  mean  value  of  the  massfraction  terms  in  Eq.  32  using  a  PDF  approach 
one  can  write 


(n„ tYmv”J)p, 


\&i (Y\ >  Yi,  ■  ■  ■ ,  Ym)dY\ ,  dY2,  ■  ■  ■ ,  dYm  (43) 


for  either  the  products  or  the  reactants.  As  the  massfractions  are  raised  to  the  power  of  the 
stoichiometric  mol  numbers  of  products  and  reactants  Eq.  43  represents  moments  of  the 
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massfraction  PDF.  Since  the  massfraction  PDF  depends  jointly  on  all  the  massfractions, 
i.e.,  the  massfractions  cannot  be  treated  as  mutually  independent, 


)W .  >2,  ■  •  • ,  r.)  ^  p»,  m  )F*(y2)  ■  ■  (44) 

this  makes  the  numerical  integration  of  even  an  assumed  Gaussian  PDF  difficult.  Therefore, 
Girimaji  [16]  proposed  the  use  of  a  multivariate  /3  distribution  defined  as, 


where  8  is  the  Dirac 


y  n.  _  r(/3i  +  pi  + . . . + fim)  \ 
~  r(Pi)r(P2)---r(pm)  l 

delta  function  and  the  variable  / 3 


r^lri2~l  ■  •  •  ^"'-1 5  (l  -  [7!  - r2  - . . .  -  Ym}) 

_  (45) 

=  f{  Ym .  Y”  Y” )  and  is  defined  as, 


Pm 


Ym 


1-X%  A 

1V>y;>  J 


(46) 


As  with  the  assumed  temperature  PDF,  this  assumed  massfraction  PDF  requires  a  knowl¬ 
edge  of  both  the  average  massfraction  value  and  a  value  for  its  variance.  Most  numerical 
models  that  consider  more  than  a  single  composition  fluid  will  have  a  species  conservation 
equation  similar  to  Eq.  10  and  so  Ym  will  be  known.  The  sum  of  the  square  of  the  fluctu¬ 
ations  in  the  species  massfractions  (sometimes  referred  to  as  the  turbulent  scalar  energy) 
is  not  typically  a  value  that  is  solved  for.  However,  development  of  variable  Schmidt  num¬ 
ber  codes  (see  Baurle  et  al.  [19])  often  solve  for  this  value  (and  the  rate  of  its  dissipation) 
as  a  consequence  of  avoiding  the  pre-specification  of  a  turbulent  Schmidt  number  to  deter¬ 
mine  the  turbulent  mass  diffusivity  (Dm)r-  In  general  the  form  of  the  massfraction  variance 
equation  can  be  expressed  as, 


dt  (P °Y )  +  dXj  (ptijOr)  —  gx.  (p  [Dm  +  CY 1  (Dm ) t\  d°J.  ) 
+2 (p(Dm) T  (#* )2  - pey  +  X 

where 


m 


(47) 


(48) 


Along  with  the  solution  of  Eq.  47,  most  variable  Schmidt  number  models  solve  for  a  similar 
equation  for  the  dissipation  rate  of  ay,  or  £y,  to  close  the  set  of  equations. 
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With  both  the  averaged  value  of  massfraction  and  its  variance  solved  for,  /3,„  is  completely 
determined  and  so  returning  to  Eq.  43  and  substituting  the  multivariate  /3  PDF  in  Eq.  45 
for  P{» ,  y2,  ■  ■  ■ ,  Ym)  yields 


(Y^JY2v^---YmVmJ) 


-  '(/ti+fe+---+fl)i  i 

r(Pi)r(ft)-r(j3U) 


n, ■  ■  ■  jZo fa V|J>2V2^  ■  •  •  7^)  (V, 


fil  ~  1  y/F 
1  22 


1  yfli. 

' ' ' 1  m 


5(i-[7t-y2 


-7m])c/7i,c/72,...,c/7m 


<7ivi./72vy-7mv»j) 


r(Pi)r(fc)-r(A») 


roo---roo(71«'72^...7m^)5(l^[71-72-...-7w])c/71,t/72,...,c/7m 


(49) 


where 


Qfoj  —  Tw./  7  /);;;  1  —  0m  1  (50) 

Examining  only  the  integral,  if  one  perfonns  the  integration  with  respect  to  Ym  first, 


7,  “'  •  •  •  7 


faw— 1 
m— 1 


Ymam8  (1  fa  —  — 


Ym )  dYm 


dY\ , . . . ,  dYm—\ 


Noting  that  the  only  time  the  Dirac  delta  function  is  non  zero  is  when  its  argument  is  zero 
one  can  define, 


7m*  =  l-(71  +  ---  +  7m_1) 

while  noting  that  S(—  x)  =  S(x)  (i.e.,  the  delta  function  is  an  even  function)  and  that  by 
virtue  of  its  sifting  property  one  can  write, 


/  g(y)S(y-c)  =  g(c) 

J  —  oo 

then  letting  g(y)  =  g(Ym )  =  7“m  and  8  (y  -  c)  =  8  ( Ym  -  Y* )  one  obtains, 


(51) 


/oo 

Yma>”8  (1  -  [7i  -  . . .  -  7m_!]  -  Ym)dYm  =  (Y*)a,n  =  ( 1  -  Yx - Ym_x)a' 

-oo 

and  thus  one  is  left  with 
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If  one  defines, 


Y**=  I-Y1-Y2 - Ym- 2 


then  the  integral  with  respect  to  7m_i  can  be  written, 


V  CX\  y^W 

^  1  ^  m— 


m—2 


YZ 


n —  1  , 
-1 


(Y**  -Ym^)a-dYm^  >  dY\ , . . . , dYr 


"m—2 


(53) 


z- 


-yOCm—l 
Im- 1 


v 


y** 

1  m 


-)amdYr 


m—  1 


Noting  that  the  sum  off  all  the  massfractions  must  equal  unity  allows  one  to  rephrase  the 
above  integral  by  observing  that, 


Y\  + 12  3 - h  Im— 2  +  Ym- 1  ! 


j  +  Ym  =  1  -  (7i  +  72  +  •  •  •  +  Ym  _2 ) 


and  thus  one  can  change  the  variable  of  integration  using  the  relations, 


and 


x  = 


Ym—  1 

y** 

1m 


Y„ 


m—  1 


V 

*  m  — 


ff«— 1 


1  Y\  Y2 - Ym 


Ym—  1  7  Ym 


(54) 


dYm_l  =Y**dx  (55) 

When  there  is  no  species  Ym_  \  then  x  =  0  while  in  the  case  where  the  entire  mixture  is 
simply  species  Ym_\  then  it  must  be  that  Ym  —  0  and  thus  x  =  1.  Therefore,  using  the  new 
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variable  of  integration  in  Eq.  54  and  applying  the  new  limits  yields  for  the  integration  with 
respect  to  Ym_  \ , 


(1  -x)a-{YZdx}  =  {YZ 


1  -|-  CCm  “b  CCm  —  1 


- x)a,ndx 


and  so  the  integral  with  respect  to  Ym-i  can  be  written, 


yam- 1  (y**  _  y  \am^V 
1m- 1  Cttj  1m- 1)  ulm 


_  /TT**\l  +  0m 

-1  -  \  /77  ) 


-l  +  0m-l  — 1 


X 


[\ -x)6m-]dx  (56) 


At  this  point  one  of  the  first  advantages  of  using  the  assumed  PDF  form  as  presented  in 
Eq.  45  presents  itself  in  that  the  remaining  integral  is  by  definition  the  beta  function  which 
itself  can  be  expressed  as  a  function  of  several  gamma  functions, 


l 


xh~l 


(1—  x)‘2  1  dx 


r(7i)r(f2) 

r(*i  +t2) 


(57) 


Combining  the  results  of  Eqs.  57  and  56  allows  the  integral  portion  of  Eq.  49  to  be  ex¬ 
pressed  as  (while  also  using  the  relation  in  Eq.  50), 


Ym\)dY],dY2,...,dYm  = 

-Ym_2)em+e'n-l~1dYu...,dYm_2 
(58)  ' 

The  integral  in  the  above  expression  has  a  form  similar  to  that  in  Eq.  52  except  that  the 
outer  two  integrations  have  been  performed.  If  one  defines, 


,/'-oo  •  •  •  ,/Too  (Y\a[  Y2ai  ■  ■  ■  Yma,n)  8  (  \  ~[Y\  —Y2  —  ... 


r(em)r(e, 


'm—l ) 


r(6m+em- 1 ) 


fit 


0,-1 


■■■Y, 


0;«- 2-1 


777  —  2 


[1-Y1-Y2---- 


Yin  -  6777  4“  0777—1  (59) 

and 

y***  \  V  V  V 

Ym  —  V—I\—I2 - Ym-3 

then  Eq.  58  can  be  expressed  as, 


r(0m)r(0,„_i)  /■“ 

r (  0777  +  0777—  1  ) 


y@m— 2 

Im- 2  J 


\*m 


Ym_2yn'-1dYh...,dYm_2 


18 


and  thus  the  integration  with  respect  to  Ym  _ 2  can  be  isolated, 


r(0m)r(0m-i) 


y  01  — 1  y0»i— 3  1 

71  ”'Im-3 


r(0m  +  0m_i)  ,/- 

and  evaluated  in  a  similar  manner  to  that  in  Eq.  53 


)  {  f  Ym—2  1  (C*  -  dYm_2  |  dY\ . . . .  .dYm_3 

(60) 


V 


/y***\  Ym~  1 
@m—2  1  /  1  m 


m— 2 


y*** 

1m 


<V*** dYm- 2 


'Em  —  ^772— 2  J 


(C*)7"'1  /  EX  1  1  - 


Ym—2 

y*** 

1m 


Ym  1 


</7 


m—2 


Noting  that 


7  E  72  3 - Ym—3  E  Ym- 2  E  Ym- 1  +  7m  —  1  —  7,***  +  7,„_2  E  7,-1  +  7m 


thus 


1  -7* 


=  Em_' 


E  7,-1  E  7 


allowing  one  to  again  change  the  variable  of  integration  by  defining, 


Ym—2 

y*** 

1  m 


Y 


m—2 


Ym~ 2  E  7,1-1  E  Yr 


and 


dYm_  2  =  Y***dx 


where  as  before  when  there  is  no  species  7„,_ 2  then  x  =  0  while  if  the  entire  mixture  is  pure 
7,-2  then  all  remaining  massfractions  must  be  zero  and  thus  x  =  I .  With  this  change  the 
integral  becomes, 


(1-x)*-1  {£**<**} 


/  y*** 

\*m 


\Ym  — 1  +  0,  ,i-2  — 1  +  1 


”m—2 


-1 


(1 


1  dx 


I  0m— 2  1 ) 


r(0„i_2)r(ym) 

r(  0„i-2  e  ) 
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where  the  definition  of  the  beta  function  in  Eq.  57  has  been  applied.  Replacing  the  holders 
ym  and  Y***  yields  for  the  integral  with  respect  to  7W_ 2, 


xr@m— 2  1  / 

Im-2  Vm 


Ym-2)7m  1  dJm- 2  =  (1  -Yi  -Y2  -•  •  •  Ym_ ^-2+e^+e, 


i)Y(6m-2)Y(dm  +  dm-i) 

Y(0m- 2  +  0m— 1  +  0m) 

(61) 


When  this  results  is  substituted  back  into  Eq.  60  the  result  becomes, 


r{Gm)r(9m-i)  r(e,„  2)170, ,,+0,,,  1) 

r(Gm+em-i)  r(em_2+em_i+0m) 


- 1 


1  y^m— 3 

Im—3 


)(!  —  7i  —  I2 - 7m-; 


(^m— l_t“0m  1) 


dYi,...,dYm_3 

(62) 


Cancelling  like  terms  it  can  now  be  noted  that  the  integrals  in  Eq.  58  and  62  have  the 
identical  fonn  but  for  the  reduction  in  the  variables  of  integration  by  one.  Therefore,  the 
process  outlined  between  Eqs.  60  and  Eq.  62  can  be  repeated  until  all  the  integrations  have 
been  performed  yielding  the  final  result  for  the  integral  in  Eq.  49  (replacing  the  holder  0m), 


JT„,  ■  ■  ■  n  (7i  yi  j'72V2j  ■  •  ■  Ymv"J)  (if1  V-f2  1---7ilm  1 )  d  (1  -  [7,  -  72  - . . .  -  7m])  J7, ,  t/72, . . .  ^7m 

r(v1./+j31)r(y;,/+fe)-r(vm>y-+ftw) 

f(fil  +&H - fAn+Vl  ,y+ V2  j-4 - hVmj) 

(63) 

With  this  integral  evaluated,  the  mean  value  of  the  massfraction  tenns  appearing  in  Eq.  32 
can  be  re-expressed  by  combining  Eq.  63  with  Eq.  49  to  obtain  the  final  result, 


/fl  V  VmJ) 
\11m1m  / 


r(/3i  +  j82  + . . .  +  pm)  T(yij  +  pi)r(v2j  +  Pi)  •  •  •  Y{vmj  +  /3,„) 
r(/3i)r(/32)  •  •  •  r(j8m)  r(/3i  +/32  h - 1-  pm  +  vij  +  v2/-  -i - 1-  vmy) 


It  is  at  this  point  that  the  second  benefit  of  using  the  assumed  form  of  the  PDF  in  Eq.  45 
becomes  apparent.  Since  the  gamma  function  F(/3 )  itself  has  the  definition, 

r(f)=  /V-'e'Vx  (65) 

Jo 


20 


one  can  make  use  of  the  fact  that  for  positive  integer  values  an  equivalent  definition  of  this 
function  is, 


r(0  =  ('-l)!  =  (l)(2)-"(<-l)  (66) 

and  thus  one  can  write, 

r(^+  1)  =  (1)(2)  •  •  •  ([*  +  1]  —  2)([*  +  1]  —  1)  =  {\){2)---{t-'\)t  =  t[{t—  1) !] 
leading  to  the  identity, 


r(r  +  i)  =  rr(0  (67) 

Repeated  use  of  this  identity  allows  the  gamma  functions  in  Eq.  64  to  be  continually  re¬ 
duced  since  vmj  is  an  integer  value, 

4"  t'mj)  =  {Pm  4"  Vm.j  1)  Y{pm  4"  V mj  1) 

' - V - - - ' 

/ - A - -s 

—  ( Pin  4“  Vmj  1)  ( Pm  4“  Ymj  2)  Y(pnl  4~  Vm,j  2) 

' - - ' 

/ - A - s 

—  {Pm  4“  Vmj  1)  {Pin  4“  YmJ  2)  {pm  4~  Ymj  3)T {Pm  4“  Ymj  2) 

Therefore, 

nPm  +  Vmj)  =  Y{Pm)U‘Zrj{Pm  +  Vmj  ~  0  (68) 

Similarly, 

r(Ift„  +  Iv,„J)  =  r(Xft„)n‘:f"v"J(Xft„  +  Xv„J-t)  (69) 

m  m  m  m  m 

Therefore,  as  these  two  terms  appear  as  a  ratio  in  Eq.  64  one  can  write, 


/n  y  Vm’j\  — 

\LLm1m  /  — 


r(2»  ft.)  (iWiCkft  +  vi  j  -  /)}  {r(j32)nj 5v(|32  +  v2J  -<)}•••  {r(ft,)n 


r(ft  )Uft)  •  •  •  r(A„)r(sm  A„ +Z,„  vmJ-k) 
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where  after  cancelling  all  the  like  terms  reduces  to  the  simple  algebraic  expression, 


/n  Y  — 

X1  Lm1m  /  — 


rG’J(Pi  +  vij 


)}{n;=rJ(j3 2  +  V2J- 


nj  ^(Pni  +  VmJ-i) 


m  VmJ  (I*  Pm  +  lm  Vmj  ~  k) 


(70) 

where  this  expression  applies  equally  to  both  the  product  and  reactant  terms.  It  is  inter¬ 
esting  to  note  that  the  solution  for  the  mean  value  of  the  terms  when  using  the  assumed 
multivariate  /3  PDF  form  in  Eq.  45  does  not  actually  require  the  solution  to  the  PDF  di¬ 
rectly.  This  is  in  contrast  to  that  for  the  temperature,  where  although  assuming  a  seemingly 
more  intuitive  Gaussian  PDF,  its  solution  requires  integration  as  shown  in  Eq.  42.  To  use 
Eq.  70  all  that  is  required  is  the  mean  value  of  the  species  massfraction  and  its  variance 
(Eq.  48),  both  of  which  can  be  found  from  conservation  equations  likely  to  already  exist  in 
a  variable  Schmidt  number  solver. 


As  a  final  consideration  it  should  be  noted  that  in  the  species  massfraction  variance  equation 
(Eq.  47)  the  tenn  Y^<bm  appears.  Re-expressing  the  species  massfraction  as  the  sum  of  its 
Favre  average  and  fluctuating  components, 


Ym(bm  —  ( Ym  ~F  ) d);;;  —  ( Ym(Om  -f  Y^ (bm) 

mm  m 

this  can  be  re-arranged  to  yield, 


'Sy' ,  Yjt[  OJ/n  —  ^  Ym  (bm  Ym  (bm 


m 


(71) 


Term  A  simply  requires  the  mean  value  for  the  chemical  source  term  which  has  already 
been  found  through  the  procedure  outlined.  For  the  second  tenn  (B)  using  the  expression 
for  the  chemical  source  term  in  Eq.  32  one  can  write, 


Ym  (bm  —  Y^  6)  p  — 

^5 1 

(72) 

where  again  the  time  averaged  massfraction  terms  have  been  replaced  with  the  expected, 
or  mean  values.  Using  the  same  assumed  form  for  the  massfraction  PDF,  the  means  ap¬ 
pearing  in  Eq.  72  are  simply  moments  of  this  PDF  and  thus  can  be  computed  in  the  same 
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fashion  as  done  for  the  means  appearing  in  Eq.  32.  In  this  case  the  moment  being  sought 
is, 


(Y^UmYmv^J)  = 


f*  oo  /»oo 


-oo  J  — oo 


{YkWmYm^)  (Ei ,  Y2, . . . ,  Ym)dYx  ,dY2,--,  dYm 

(73) 

where  after  substituting  the  assumed  form  for  the  PDF  in  Eq.  45  yields, 


/V,  Y  V|  ,/'  Y-,  V2J  ■  ■  ■  Y  V>nj\  —  n/jl+fe  +  --+fi 

vs1 1  12  lm  >  -  m)nih)-npm 


roo---Sroo(Y^Y2^---Yma'n)8(l-[Yl-Y2-...-Ym])dY1,dY2,...,dYm 


(74) 


where  ccm  has  the  same  definition  as  in  Eq.  50.  Perfonning  the  same  steps  as  outlined  from 
Eq.  49  for  all  the  variables  of  integration  except  for  dY^  will  yield  the  same  result  as  shown 
in  Eq.  62  leaving  for  the  final  integration  step, 


r(e1)r(02)---r(ow_1)  r°°  ,  0m-i 
r{ol  +  e2  +  ~-  +  em-l )J—\  m  m 


(1  -Ym)^+e2+-+6m-^dYm 


Defining 


^m  —  0i+02-\ - b  Qm - 1  and  x  =  Ym  and  thus  dx  —  dYm  (75) 

allows  the  integral  to  be  rephrased  as, 


r(e1)r(e2)---r(efW_1) 
r(0i  +  o2-\ - f  o,n- 1) 


xGm(l-x)^m  ^dx 


where  as  before,  if  the  mixture  has  no  Ym  then  x  =  0  while  for  a  mixture  of  pure  Ym,  x  —  1 . 
Using  the  result  for  the  beta  function  (Eq.  57)  where  in  this  case  t\  —  1  =  Qm  allows  one  to 
replace  the  integral  to  obtain, 


r(0i)r(o2) •  •  ■  r(om_i)  r(em  +  i)r(^) 

r(@i  +  02  4 - b  Om-l)  r (6m  +  1  +  %rn) 


Noting  that  T(E,m)  —  r(@i  +  02  H - b  6m-i)  these  terms  cancel  leaving, 


r(0i)r(o2)  •  •  •  r(em_i)r(fl„,  + 1) 
r(0/7!  +  1  +  ^ m ) 


(76) 
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From  the  identity  expressed  in  Eq.  67  one  can  write 


r  (dm  + 1)  —  o,n  r(0,„) 


and 


r(<^7 n  +  Qm  +  1)  —  r(01  +  02  H - 6  0/72-1  +  0/72  + 1) 


—  (01  +  02  H - 6  0/72-1  +  0m)r(0i  +  02  H - b  0/72-1  +  0m) 

which  allows  Eq.  76  to  be  written, 


0m{r(01)r(02)---r(0m 1)r(0m)} 


(77) 


(01  +  02  H - b  0„,)r(0i  +  02  H - b  0/72 ) 


This  has  the  exact  same  form  as  the  result  expressed  in  Eq.  63  but  for  the  additional  pre- 
multiplying  term, 


iym,j  “b  Pm) 


(78) 


(01  +  02  H - b  0/22- 1  +  0rn)  (Pi  +  p2  H - b  Pm  +  V]  j  +  V2  j  H - b  Vmj) 


and  thus  application  of  the  results  expressed  by  Eqs.  68  and  69  will  yield  a  final  result 
almost  identical  to  that  found  in  Eq.  70, 


(79) 


where  again  this  applies  equally  to  the  product  and  reactant  terms  appearing  in  Eq.  72.  As 
with  Eq.  70  the  solution  to  the  assumed  multivariate  P  PDF  is  not  required  to  evaluate  the 
mean  as  expressed  in  Eq.  79.  As  before,  if  one  can  calculate  p  =  f(Ym ,  oy)  (Eq.  46)  then 
along  with  the  stoichiometric  mol  numbers  of  the  reactant  and  product  species  involved 
in  a  given  reaction  j,  one  need  only  solve  an  algebraic  expression  to  obtain  the  time  aver¬ 
aged  chemical  source  terms  appearing  in  the  species  continuity  (Eq.  10)  and  massfraction 
variance  (Eq.  47)  equations. 
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3  Recommendations 


Due  to  its  resulting  computational  simplicity,  the  use  of  an  assumed  multivariate  /3  PDF  has 
been  used  in  a  wide  variety  of  flowfields  outside  the  original  scope  of  its  original  proposed 
use  for  turbulent  mixing.  By  assuming  statistical  independence  between  the  temperature 
and  the  massfraction  fluctuations,  the  assumed  massfraction  PDF  matches  the  behaviour 
of  the  real  PDF  in  regions  where  turbulent  mixing  alone  occurs.  It  can  be  shown  that  the 
behaviour  of  this  assumed  PDF  in  the  limit  approaches  that  of  a  standardized  Gaussian 
PDF  (which  matches  observed  results  from  Direct  Numerical  Simulations  for  two  scalar 
mixing  over  all  stages  of  the  mixing  process).  However  the  highly  dissipative  nature  of 
the  tenn  F;"  d)m  found  in  the  conservation  equation  for  the  dissipation  of  the  sum  of  the 
massfraction  variance  causes  discrepancies  between  the  model  and  experiment  for  regions 
where  chemical  reactions  are  occurring  (mean  flow  properties  are  reasonably  consistent 
with  experiment).  This  has  led  recent  researchers  to  simply  model  these  terms  to  improve 
comparison  with  experiment  (Keistler,  Hassan,  and  Xiao  [29], [30],  Keistler  and  Hassan 
[31]). 

Evolution  PDFs,  although  capable  of  potentially  overcoming  some  of  the  issues  with  the 
assumed  PDF  approach  by  incorporating  more  variables  within  a  single  joint  PDF,  still 
require  computational  resources  that  make  their  use  prohibitive.  However,  their  strong  the¬ 
oretical  basis  make  them  an  attractive  alternative  and  thus  it  is  recommended  a  comprehen¬ 
sive  review  of  evolution  PDF  methods  and  their  recent  application  to  high  speed  reacting 
flows  be  conducted.  Within  this  review,  alternatives  to  the  assumed  multivariate  /3  PDF 
should  be  examined  in  an  effort  to  establish  if  progress  has  been  made  to  reduce  the  dissi¬ 
pative  effect  of  certain  terms  that  appear  within  the  variance  dissipation  equation. 
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List  of  Acronyms  and  Symbols 


Acronyms 

PDF  Probability  Density  Function 

DRDC  Defence  Research  and  Development  Canada 


Greek  Symbols 

8  Dirac  delta  function 

£y  Dissipation  of  massfraction  variance 

r  Gamma  function 

K  Laminar  thermal  conductivity 

p  Laminar  co-efficient  of  viscosity 

vm  Species  mass  diffusivity,  stoichiometric  mole  number 

p  Density 

a  Standard  deviation 

ay  Sum  of  mass  fraction  variances 

co  Specific  dissipation  of  k 

(bm  Rate  of  production  or  consumption  of  species  m 

£  Enstrophy 


Roman  Symbols 

Cp  Specific  heat  at  constant  pressure 

Cv  Specific  heat  at  constant  volume 

Dm  Species  mass  diffusion  co-efficient 

e  Specific  energy 

E  Expectation 

Ea  Activation  energy 

h  Specific  enthalpy 

k  Specific  turbulent  kinetic  energy 

kb  Backward  reaction  rate 
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kf  Forward  reaction  rate 

K  Reaction  rate 

Kc  Equilibrium  constant  based  on  concentrations 

.JYm  Species  molecular  weight 

N  Number  of  experiments 

Probability  density  function  of  (p 
£%  Universal  gas  constant 

Sc  Schmidt  number 

t  Time 

T  Temperature,  time 

u  velocity 

Y  Mass  fraction 


Subscripts 

m  Species 

T  Turbulent 


Superscripts 

Favre  fluctuating  component 
p  Product 

r  Reactant 
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